suppressPackageStartupMessages({
require(tidyr)
require(lme4) # for running models
require(sjPlot) # for visualizing model interactions
require(sjmisc)
require(sjlabelled)
require(tidyverse)
require(viridis)
require(lmerTest)
})
Population_climateNA <- read.csv("C:/Dropbox/UVM/Research/Code/ClimateData/ClimateNA/Family/PCA_sprucePops.txt",
header=T, sep="\t")
colnames(Population_climateNA)[colnames(Population_climateNA)=="PC1_sel"] <- "PC1"
Family_climateNA <- read.csv("C:/Dropbox/UVM/Research/Code/ClimateData/ClimateNA/Family/PCA_spruceFams.txt",
header=T, sep="\t")
colnames(Family_climateNA)[colnames(Family_climateNA)=="PC1_sel"] <- "PC1"
budset_2019 <- read.csv("./trait_data/BudSet_2019.csv", stringsAsFactors = T)
budset_2020 <- read.csv("./trait_data/BudSet_2020.csv", stringsAsFactors = T)
budbreak_2020 <- read.csv("./trait_data/BudBreak_2020_cGDD.csv", stringsAsFactors = T)
growth_2019 <- read.csv("./trait_data/Growth_2019.csv", stringsAsFactors = T)
growth_2020 <- read.csv("./trait_data/Growth_2020.csv", stringsAsFactors = T)
Plasticity <- read.csv("./trait_data/Plasticity_TD.csv", stringsAsFactors = T)
# Garden
budset_2019$Garden <- factor(budset_2019$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
budset_2020$Garden <- factor(budset_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
budbreak_2020$Garden <- factor(budbreak_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
growth_2019$Garden <- factor(growth_2019$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
growth_2020$Garden <- factor(growth_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
# Region
budset_2019$Region <- factor(budset_2019$Region,levels = c("Core", "Margin", "Edge"))
budset_2020$Region <- factor(budset_2020$Region,levels = c("Core", "Margin", "Edge"))
budbreak_2020$Region <- factor(budbreak_2020$Region,levels = c("Core", "Margin", "Edge"))
growth_2019$Region <- factor(growth_2019$Region,levels = c("Core", "Margin", "Edge"))
growth_2020$Region <- factor(growth_2020$Region,levels = c("Core", "Margin", "Edge"))
## working model for growth - ****
Growth19_mod1 <- lmer(data = growth_2019, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E
## boundary (singular) fit: see ?isSingular
Growth19_mod1a <- lmer(data = growth_2019, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
Growth19_mod1b <- lmer(data = growth_2019, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(Growth19_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 75.68 37.84 2 12.32 1.0818 0.369101
## Region 2991.90 1495.95 2 64.21 42.7664 1.562e-12 ***
## Garden:Region 532.33 133.08 4 722.10 3.8046 0.004536 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(Growth19_mod1, Growth19_mod1a)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## Growth19_mod1a: Family)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## Growth19_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1a 13 30358 30442 -15166 30332
## Growth19_mod1 18 30351 30467 -15158 30315 16.499 5 0.005555 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(Growth19_mod1a, Growth19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## Growth19_mod1a: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1b 12 30469 30547 -15223 30445
## Growth19_mod1a 13 30358 30442 -15166 30332 113.39 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(Growth19_mod1, Growth19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## Growth19_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1b 12 30469 30547 -15223 30445
## Growth19_mod1 18 30351 30467 -15158 30315 129.89 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(Growth19_mod1, Growth19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## Growth19_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1b 12 30469 30547 -15223 30445
## Growth19_mod1 18 30351 30467 -15158 30315 129.89 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT
anova(Growth19_mod1,Growth19_mod1a,Growth19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## Growth19_mod1a: Family)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## Growth19_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1b 12 30469 30547 -15223 30445
## Growth19_mod1a 13 30358 30442 -15166 30332 113.392 1 < 2.2e-16 ***
## Growth19_mod1 18 30351 30467 -15158 30315 16.499 5 0.005555 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(Growth19_mod1,Growth19_mod1a,Growth19_mod1b)
## df AIC
## Growth19_mod1 18 30344.87
## Growth19_mod1a 13 30351.38
## Growth19_mod1b 12 30462.51
# effect size estimates for growth trait
plot_model(Growth19_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
# diagnostic tests for the model
# sjPlot::plot_model(Growth19_mod1, type = "diag")
tab_model(Growth19_mod1)
| Growth | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 9.53 | 7.86 – 11.20 | <0.001 |
| Garden [Maryland] | 0.03 | -2.24 – 2.31 | 0.976 |
| Garden [North_Carolina] | -0.41 | -2.69 – 1.86 | 0.721 |
| Region [Margin] | 4.74 | 3.49 – 6.00 | <0.001 |
| Region [Edge] | 3.02 | 2.05 – 3.98 | <0.001 |
|
Garden [Maryland] * Region [Margin] |
0.41 | -0.79 – 1.61 | 0.506 |
|
Garden [North_Carolina] * Region [Margin] |
-1.73 | -2.95 – -0.51 | 0.005 |
|
Garden [Maryland] * Region [Edge] |
0.29 | -0.65 – 1.23 | 0.548 |
|
Garden [North_Carolina] * Region [Edge] |
-0.96 | -1.92 – -0.00 | 0.050 |
| Random Effects | |||
| σ2 | 34.98 | ||
| τ00 Family | 3.55 | ||
| τ00 Population | 1.03 | ||
| τ00 mBed | 3.13 | ||
| τ11 Family.GardenMaryland | 0.89 | ||
| τ11 Family.GardenNorth_Carolina | 0.47 | ||
| ρ01 Family.GardenMaryland | 0.35 | ||
| ρ01 Family.GardenNorth_Carolina | -0.95 | ||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4681 | ||
| Marginal R2 / Conditional R2 | 0.089 / NA | ||
## working model for growth - ****
Growth_mod1 <- lmer(data = growth_2020, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -2.3e+01
Growth_mod1a <- lmer(data = growth_2020, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
Growth_mod1b <- lmer(data = growth_2020, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(Growth_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 1602.55 801.27 2 13.77 17.6768 0.0001573 ***
## Region 2830.07 1415.03 2 67.79 31.2169 2.454e-10 ***
## Garden:Region 815.96 203.99 4 368.64 4.5002 0.0014624 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(Growth_mod1, Growth_mod1a)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## Growth_mod1a: Family)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## Growth_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1a 13 30161 30244 -15067 30135
## Growth_mod1 18 30164 30280 -15064 30128 6.1533 5 0.2916
# effect of genotype
anova(Growth_mod1a, Growth_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## Growth_mod1a: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1b 12 30216 30293 -15096 30192
## Growth_mod1a 13 30161 30244 -15067 30135 57.164 1 4.009e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(Growth_mod1, Growth_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## Growth_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1b 12 30216 30293 -15096 30192
## Growth_mod1 18 30164 30280 -15064 30128 63.317 6 9.51e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(Growth_mod1, Growth_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## Growth_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1b 12 30216 30293 -15096 30192
## Growth_mod1 18 30164 30280 -15064 30128 63.317 6 9.51e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT
anova(Growth_mod1,Growth_mod1a,Growth_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## Growth_mod1a: Family)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## Growth_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1b 12 30216 30293 -15096 30192
## Growth_mod1a 13 30161 30244 -15067 30135 57.1641 1 4.009e-14 ***
## Growth_mod1 18 30164 30280 -15064 30128 6.1533 5 0.2916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(Growth_mod1,Growth_mod1a,Growth_mod1b)
## df AIC
## Growth_mod1 18 30157.95
## Growth_mod1a 13 30154.50
## Growth_mod1b 12 30209.53
# effect size estimates for growth trait
plot_model(Growth_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
# diagnostic tests for the model
# sjPlot::plot_model(Growth19_mod1, type = "diag")
tab_model(Growth_mod1)
| Growth | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 16.04 | 14.91 – 17.17 | <0.001 |
| Garden [Maryland] | 0.87 | -0.50 – 2.25 | 0.212 |
| Garden [North_Carolina] | -2.94 | -4.33 – -1.56 | <0.001 |
| Region [Margin] | 3.76 | 2.16 – 5.37 | <0.001 |
| Region [Edge] | -0.16 | -1.39 – 1.07 | 0.796 |
|
Garden [Maryland] * Region [Margin] |
3.08 | 1.56 – 4.59 | <0.001 |
|
Garden [North_Carolina] * Region [Margin] |
2.69 | 1.12 – 4.26 | 0.001 |
|
Garden [Maryland] * Region [Edge] |
0.49 | -0.71 – 1.69 | 0.420 |
|
Garden [North_Carolina] * Region [Edge] |
0.67 | -0.57 – 1.92 | 0.289 |
| Random Effects | |||
| σ2 | 45.33 | ||
| τ00 Family | 0.00 | ||
| τ00 Population | 3.10 | ||
| τ00 mBed | 0.86 | ||
| τ11 Family.GardenMaryland | 4.07 | ||
| τ11 Family.GardenNorth_Carolina | 6.30 | ||
| ρ01 Family.GardenMaryland | |||
| ρ01 Family.GardenNorth_Carolina | |||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4475 | ||
| Marginal R2 / Conditional R2 | 0.136 / NA | ||
budset_2019 <- budset_2019%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
## working model for growth - ****
budset19_mod1 <- lmer(data = budset_2019, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## Warning: Model failed to converge with 1 negative eigenvalue: -3.2e+01
budset19_mod1a <- lmer(data = budset_2019, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
budset19_mod1b <- lmer(data = budset_2019, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(budset19_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 39470 19735.0 2 13.60 31.1200 8.408e-06 ***
## Region 39318 19658.9 2 69.07 30.9999 2.465e-10 ***
## Garden:Region 12902 3225.6 4 604.80 5.0865 0.0004895 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(budset19_mod1, budset19_mod1a)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budset19_mod1a: Family)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budset19_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1a 13 44060 44144 -22017 44034
## budset19_mod1 18 44019 44135 -21991 43983 50.972 5 8.765e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(budset19_mod1a, budset19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budset19_mod1a: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1b 12 44135 44212 -22055 44111
## budset19_mod1a 13 44060 44144 -22017 44034 77.107 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(budset19_mod1, budset19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budset19_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1b 12 44135 44212 -22055 44111
## budset19_mod1 18 44019 44135 -21991 43983 128.08 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(budset19_mod1, budset19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budset19_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1b 12 44135 44212 -22055 44111
## budset19_mod1 18 44019 44135 -21991 43983 128.08 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT
anova(budset19_mod1,budset19_mod1a,budset19_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budset19_mod1a: Family)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budset19_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1b 12 44135 44212 -22055 44111
## budset19_mod1a 13 44060 44144 -22017 44034 77.107 1 < 2.2e-16 ***
## budset19_mod1 18 44019 44135 -21991 43983 50.972 5 8.765e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(budset19_mod1,budset19_mod1a,budset19_mod1b)
## df AIC
## budset19_mod1 18 43989.61
## budset19_mod1a 13 44029.72
## budset19_mod1b 12 44104.62
# effect size estimates for budset trait
plot_model(budset19_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
tab_model(budset19_mod1)
| BudSet | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 260.61 | 255.82 – 265.40 | <0.001 |
| Garden [Maryland] | -8.61 | -14.69 – -2.52 | 0.006 |
| Garden [North_Carolina] | -20.62 | -26.83 – -14.40 | <0.001 |
| Region [Margin] | 8.89 | 3.12 – 14.66 | 0.003 |
| Region [Edge] | 16.94 | 12.49 – 21.39 | <0.001 |
|
Garden [Maryland] * Region [Margin] |
-4.30 | -9.64 – 1.05 | 0.115 |
|
Garden [North_Carolina] * Region [Margin] |
-8.88 | -15.02 – -2.74 | 0.005 |
|
Garden [Maryland] * Region [Edge] |
-7.34 | -11.55 – -3.13 | 0.001 |
|
Garden [North_Carolina] * Region [Edge] |
-1.85 | -6.66 – 2.96 | 0.450 |
| Random Effects | |||
| σ2 | 634.16 | ||
| τ00 Family | 103.15 | ||
| τ00 Population | 20.11 | ||
| τ00 mBed | 19.75 | ||
| τ11 Family.GardenMaryland | 42.42 | ||
| τ11 Family.GardenNorth_Carolina | 100.27 | ||
| ρ01 Family.GardenMaryland | -1.00 | ||
| ρ01 Family.GardenNorth_Carolina | -0.39 | ||
| ICC | 0.16 | ||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4685 | ||
| Marginal R2 / Conditional R2 | 0.143 / 0.277 | ||
budset_2020 <- budset_2020%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
## working model for growth - ****
budset20_mod1 <- lmer(data = budset_2020, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E
## boundary (singular) fit: see ?isSingular
budset20_mod1a <- lmer(data = budset_2020, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
budset20_mod1b <- lmer(data = budset_2020, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(budset20_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 66924 33462 2 13.29 107.4274 6.252e-09 ***
## Region 6295 3148 2 62.25 10.1050 0.0001583 ***
## Garden:Region 3817 954 4 518.61 3.0635 0.0163764 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(budset20_mod1, budset20_mod1a)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budset20_mod1a: Family)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budset20_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1a 13 37141 37224 -18558 37115
## budset20_mod1 18 37116 37231 -18540 37080 35.178 5 1.386e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(budset20_mod1a, budset20_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budset20_mod1a: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1b 12 37185 37261 -18580 37161
## budset20_mod1a 13 37141 37224 -18558 37115 45.36 1 1.64e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(budset20_mod1, budset20_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budset20_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1b 12 37185 37261 -18580 37161
## budset20_mod1 18 37116 37231 -18540 37080 80.538 6 2.766e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(budset20_mod1, budset20_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budset20_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1b 12 37185 37261 -18580 37161
## budset20_mod1 18 37116 37231 -18540 37080 80.538 6 2.766e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT
anova(budset20_mod1,budset20_mod1a,budset20_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budset20_mod1a: Family)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budset20_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1b 12 37185 37261 -18580 37161
## budset20_mod1a 13 37141 37224 -18558 37115 45.360 1 1.640e-11 ***
## budset20_mod1 18 37116 37231 -18540 37080 35.178 5 1.386e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(budset20_mod1,budset20_mod1a,budset20_mod1b)
## df AIC
## budset20_mod1 18 37092.72
## budset20_mod1a 13 37118.40
## budset20_mod1b 12 37161.59
# effect size estimates for budset trait
plot_model(budset20_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
tab_model(budset20_mod1)
| BudSet | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 235.37 | 232.04 – 238.70 | <0.001 |
| Garden [Maryland] | -25.29 | -29.69 – -20.89 | <0.001 |
| Garden [North_Carolina] | -28.81 | -33.25 – -24.37 | <0.001 |
| Region [Margin] | 6.71 | 2.89 – 10.53 | 0.001 |
| Region [Edge] | 2.02 | -0.94 – 4.97 | 0.181 |
|
Garden [Maryland] * Region [Margin] |
0.02 | -4.04 – 4.08 | 0.991 |
|
Garden [North_Carolina] * Region [Margin] |
-2.36 | -6.50 – 1.79 | 0.266 |
|
Garden [Maryland] * Region [Edge] |
-2.15 | -5.34 – 1.04 | 0.187 |
|
Garden [North_Carolina] * Region [Edge] |
2.00 | -1.33 – 5.33 | 0.240 |
| Random Effects | |||
| σ2 | 311.48 | ||
| τ00 Family | 57.16 | ||
| τ00 Population | 4.93 | ||
| τ00 mBed | 10.02 | ||
| τ11 Family.GardenMaryland | 37.27 | ||
| τ11 Family.GardenNorth_Carolina | 38.27 | ||
| ρ01 Family.GardenMaryland | -0.97 | ||
| ρ01 Family.GardenNorth_Carolina | -0.78 | ||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4282 | ||
| Marginal R2 / Conditional R2 | 0.361 / NA | ||
budbreak_2020 <- budbreak_2020%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
## working model for growth - ****
budbreak_mod1 <- lmer(data = budbreak_2020, na.action = na.omit,
cGDD ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 2 negative eigenvalues: -6.2e-02 -6.7e+01
budbreak_mod1a <- lmer(data = budbreak_2020, na.action = na.omit,
cGDD ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
budbreak_mod1b <- lmer(data = budbreak_2020, na.action = na.omit,
cGDD ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(budbreak_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 2521553 1260777 2 17.95 1104.830 < 2.2e-16 ***
## Region 122334 61167 2 90.15 53.601 4.592e-16 ***
## Garden:Region 45697 11424 4 370.29 10.011 1.056e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(budbreak_mod1, budbreak_mod1a)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budbreak_mod1a: Family)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budbreak_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1a 13 47668 47752 -23821 47642
## budbreak_mod1 18 47524 47641 -23744 47488 153.5 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(budbreak_mod1a, budbreak_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budbreak_mod1a: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1b 12 47825 47903 -23901 47801
## budbreak_mod1a 13 47668 47752 -23821 47642 159.29 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(budbreak_mod1, budbreak_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budbreak_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1b 12 47825 47903 -23901 47801
## budbreak_mod1 18 47524 47641 -23744 47488 312.79 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(budbreak_mod1, budbreak_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budbreak_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1b 12 47825 47903 -23901 47801
## budbreak_mod1 18 47524 47641 -23744 47488 312.79 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT
anova(budbreak_mod1,budbreak_mod1a,budbreak_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budbreak_mod1a: Family)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budbreak_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1b 12 47825 47903 -23901 47801
## budbreak_mod1a 13 47668 47752 -23821 47642 159.29 1 < 2.2e-16 ***
## budbreak_mod1 18 47524 47641 -23744 47488 153.50 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(budbreak_mod1,budbreak_mod1a,budbreak_mod1b)
## df AIC
## budbreak_mod1 18 47491.26
## budbreak_mod1a 13 47635.36
## budbreak_mod1b 12 47792.33
# effect size estimates for budbreak trait
plot_model(budbreak_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
tab_model(budbreak_mod1)
| cGDD | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 562.01 | 557.99 – 566.04 | <0.001 |
| Garden [Maryland] | 90.76 | 85.39 – 96.13 | <0.001 |
| Garden [North_Carolina] | 120.18 | 113.90 – 126.45 | <0.001 |
| Region [Margin] | 7.47 | 1.55 – 13.39 | 0.013 |
| Region [Edge] | 12.03 | 7.48 – 16.59 | <0.001 |
|
Garden [Maryland] * Region [Margin] |
-1.43 | -8.77 – 5.90 | 0.702 |
|
Garden [North_Carolina] * Region [Margin] |
-6.34 | -16.28 – 3.60 | 0.211 |
|
Garden [Maryland] * Region [Edge] |
10.39 | 4.63 – 16.15 | <0.001 |
|
Garden [North_Carolina] * Region [Edge] |
21.99 | 14.18 – 29.80 | <0.001 |
| Random Effects | |||
| σ2 | 1141.15 | ||
| τ00 Family | 0.00 | ||
| τ00 Population | 24.59 | ||
| τ00 mBed | 10.48 | ||
| τ11 Family.GardenMaryland | 102.53 | ||
| τ11 Family.GardenNorth_Carolina | 560.79 | ||
| ρ01 Family.GardenMaryland | |||
| ρ01 Family.GardenNorth_Carolina | |||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4757 | ||
| Marginal R2 / Conditional R2 | 0.725 / NA | ||
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
Growth19_mod_nMb <- lmer(data = growth_2019,
Growth ~ Garden * Region +
(1|Population) +
(Garden|Family))
## boundary (singular) fit: see ?isSingular
heightData19 <- growth_2019
heightData19$pred <- predict(Growth19_mod_nMb)
# merge climateNA data
heightData19 <- merge(x=heightData19,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
G19 <- ggplot(heightData19,aes(Garden,Growth))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
G19 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("Garden site") + ylab("Growth in 2019 (in cm)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
Region plot
# Rxn - population
growth19_pop_rxn <- G19 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("Vertical height growth (in cm)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(heightData19,
aes(x=Garden, y=Growth, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Growth in 2019 (in cm)") + theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
Growth20_mod_nMb <- lmer(data = growth_2020,
Growth ~ Garden * Region +
(1|Population) +
(Garden|Family))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.2e+02
heightData20 <- growth_2020
heightData20$pred <- predict(Growth20_mod_nMb)
# merge climateNA data
heightData20 <- merge(x=heightData20,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
G20 <- ggplot(heightData20,aes(Garden,Growth))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
G20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("Garden site") + ylab("Growth in 2020 (in cm)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
Region plot
# Rxn - population
growth20_pop_rxn <- G20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("Vertical height growth (in cm)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(heightData20,
aes(x=Garden, y=Growth, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Growth in 2020 (in cm)") + theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budset19_mod_nMb <- lmer(data = budset_2019,
BudSet ~ Garden * Region +
(1|Population) +
(Garden|Family))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -8.5e+00
budset_2019_2 <- budset_2019
budset_2019_2$pred <- predict(budset19_mod_nMb)
# merge climateNA data
budset_2019_2 <- merge(x=budset_2019_2,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BS19 <- ggplot(budset_2019_2,aes(Garden,BudSet))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
BS19 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("Garden site") + ylab("Bud set in 2019 (DOY)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
Region plot
# Rxn - population
budset19_pop_rxn <- BS19 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("Bud set in 2019 (DOY)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(budset_2019_2,
aes(x=Garden, y=BudSet, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Bud set in 2019 (DOY)") + theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budset20_mod_nMb <- lmer(data = budset_2020,
BudSet ~ Garden * Region +
(1|Population) +
(Garden|Family))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -6.8e+01
budset_2020_2 <- budset_2020
budset_2020_2$pred <- predict(budset20_mod_nMb)
# merge climateNA data
budset_2020_2 <- merge(x=budset_2020_2,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BS20 <- ggplot(budset_2020_2,aes(Garden,BudSet))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
BS20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("Garden site") + ylab("Bud set in 2020 (DOY)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
Region plot
# Rxn - population
budset20_pop_rxn <- BS20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("Bud set in 2020 (DOY)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(budset_2020_2,
aes(x=Garden, y=BudSet, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Bud set in 2020 (DOY)") + theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budbreak_mod_nMb <- lmer(data = budbreak_2020,
cGDD ~ Garden * Region +
(1|Population) +
(Garden|Family))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.3e+02
budbreak <- budbreak_2020[!is.na(budbreak_2020$cGDD),]
budbreak$pred <- predict(budbreak_mod_nMb)
# merge climateNA data
budbreak <- merge(x=budbreak,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BB20 <- ggplot(budbreak,aes(Garden,cGDD))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
BB20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("Garden site") + ylab("Bud break in 2020 (cGDD)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
Region plot
# Rxn - population
budbreak_pop_rxn <- BB20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("cGDD") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
ggplot(budbreak,
aes(x=Garden, y=cGDD, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Bud break in 2020 (cGDD)") + theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
budbreak_cCDD <- read.csv("C:\\Dropbox\\UVM\\Research\\Code/CommonGardenTraits/2020/BudBreak_2020_cCDD.csv", header=T,
stringsAsFactors = T)
# reordering factors
budbreak_cCDD$Region <- factor(budbreak_cCDD$Region,
levels = c("Core", "Margin", "Edge"))
budbreak_cCDD$Garden <- factor(budbreak_cCDD$Garden,
levels = c("Vermont", "Maryland", "North_Carolina"))
budbreak_cCDD <- budbreak_cCDD%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
## working model for cCDD - ****
budbreak_cCDD_mod1 <- lmer(data = budbreak_cCDD,
cCDD ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family), na.action = na.omit) # G x E
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 2 negative eigenvalues: -8.8e-02 -9.2e-01
budbreak_cCDD_mod1a <- lmer(data = budbreak_cCDD,
cCDD ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family), na.action = na.omit)
budbreak_cCDD_mod1b <- lmer(data = budbreak_cCDD,
cCDD ~ Garden * Region +
(1|Population) + (1|mBed), na.action = na.omit)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.09547 (tol = 0.002, component 1)
# model 1 anova
anova(budbreak_cCDD_mod1)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 125141 62570 2 21.26 83756.3158 < 2.2e-16 ***
## Region 34 17 2 400.59 22.5735 5.136e-10 ***
## Garden:Region 29 7 4 485.70 9.5912 1.801e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# impact of GxE
anova(budbreak_cCDD_mod1, budbreak_cCDD_mod1a)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1a: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budbreak_cCDD_mod1a: Family)
## budbreak_cCDD_mod1: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budbreak_cCDD_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_cCDD_mod1a 13 13162 13246 -6568.0 13136
## budbreak_cCDD_mod1 18 12728 12845 -6346.3 12692 443.39 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# effect of genotype
anova(budbreak_cCDD_mod1a, budbreak_cCDD_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1b: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_cCDD_mod1a: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budbreak_cCDD_mod1a: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_cCDD_mod1b 12 13230 13308 -6603.1 13206
## budbreak_cCDD_mod1a 13 13162 13246 -6568.0 13136 70.354 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 1 vs 3
anova(budbreak_cCDD_mod1, budbreak_cCDD_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1b: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_cCDD_mod1: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budbreak_cCDD_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_cCDD_mod1b 12 13230 13308 -6603.1 13206
## budbreak_cCDD_mod1 18 12728 12845 -6346.3 12692 513.74 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(budbreak_cCDD_mod1, budbreak_cCDD_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1b: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_cCDD_mod1: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budbreak_cCDD_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_cCDD_mod1b 12 13230 13308 -6603.1 13206
## budbreak_cCDD_mod1 18 12728 12845 -6346.3 12692 513.74 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT
anova(budbreak_cCDD_mod1,budbreak_cCDD_mod1a,budbreak_cCDD_mod1b)
## refitting model(s) with ML (instead of REML)
## Data: budbreak_cCDD
## Models:
## budbreak_cCDD_mod1b: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_cCDD_mod1a: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 |
## budbreak_cCDD_mod1a: Family)
## budbreak_cCDD_mod1: cCDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden |
## budbreak_cCDD_mod1: Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_cCDD_mod1b 12 13230 13308 -6603.1 13206
## budbreak_cCDD_mod1a 13 13162 13246 -6568.0 13136 70.354 1 < 2.2e-16 ***
## budbreak_cCDD_mod1 18 12728 12845 -6346.3 12692 443.389 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AIC
AIC(budbreak_cCDD_mod1,budbreak_cCDD_mod1a,budbreak_cCDD_mod1b)
## df AIC
## budbreak_cCDD_mod1 18 12762.19
## budbreak_cCDD_mod1a 13 13196.03
## budbreak_cCDD_mod1b 12 13264.19
# effect size estimates for budbreak_cCDD trait
plot_model(budbreak_cCDD_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
# reaction norms for interactions
sjPlot::plot_model(budbreak_cCDD_mod1, type = "int")
sjPlot::plot_model(budbreak_cCDD_mod1, type = "re")
## [[1]]
##
## [[2]]
##
## [[3]]
##### reaction norms based on the model - Family
budbreak_cCDD3 <- budbreak_cCDD[!is.na(budbreak_cCDD$cCDD),]
# budbreak_cCDD2$pred <- predict(budbreak_cCDD_mod1, re.form=NULL) # all random effects are included
budbreak_cCDD3$pred <- predict(budbreak_cCDD_mod1,re.form=NA) ## population level
budbreak_cCDD3$pred1 <- predict(budbreak_cCDD_mod1) ## individual level
g0 <- ggplot(budbreak_cCDD3,aes(Garden,cCDD, label=Family))+
geom_point()+
geom_line(aes(group=Family)) +
facet_grid(.~Region)
# zero_margin
g0 + geom_line(colour="blue",aes(y=pred1,group=Family)) +
geom_line(colour="red",aes(y=pred,group=Family))
##### reaction norms based on the model - Population
g0 <- ggplot(budbreak_cCDD3,aes(Garden,cCDD))+
geom_point()+
geom_line(aes(group=Population)) +
facet_grid(.~Region)
# zero_margin
g0 + geom_line(colour="blue",aes(y=pred1,group=Population)) +
geom_line(colour="red",aes(y=pred,group=Population))
# reaction norms based on the model - Family
budbreak_cCDD3 <- budbreak_cCDD[!is.na(budbreak_cCDD$cCDD),]
# budbreak_cCDD2$pred <- predict(budbreak_cCDD_mod1, re.form=NULL) # all random effects are included
budbreak_cCDD3$pred <- predict(budbreak_cCDD_mod1,re.form=NA) ## population level/ no random effects included
budbreak_cCDD3$pred1 <- predict(budbreak_cCDD_mod1) ## individual level
G0 <- ggplot(budbreak_cCDD3,aes(Garden,cCDD))+
# geom_point()+
# geom_line(aes(group=Family)) +
facet_grid(.~Region) #+ zero_margin
G0 + geom_line(colour="red",aes(y=pred,group=Family)) +
xlab("Garden site") + ylab("Vertical height growth (in cm)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
budbreak_cCDD 2020 rxn - paper
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budbreak_cCDD_mod_nMb <- lmer(data = budbreak_cCDD,
cCDD ~ Garden * Region +
(1|Population) +
(Garden|Family), na.action = na.omit)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 2 negative eigenvalues: -1.0e+00 -1.1e+00
budbreak_cCDD6 <- budbreak_cCDD3
budbreak_cCDD6$pred <- predict(budbreak_cCDD_mod_nMb)
# merge climateNA data
budbreak_cCDD6 <- merge(x=budbreak_cCDD6,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
G0 <- ggplot(budbreak_cCDD6,aes(Garden,cCDD))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
G0 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("Garden site") + ylab("Bud break in 2020 (cCDD)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
# Rxn - population
budbreak_cCDD_pop_rxn <- G0 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("cCDD") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
plotly::ggplotly(budbreak_cCDD_pop_rxn)
budbreak_cCDD_pop_rxn
## Reaction norm for geo genetic groups/regions
# growth reaction norm
ggplot(budbreak_cCDD3, aes(x = Garden, y = cCDD, color = Region)) +
stat_summary(aes(group = Region), fun = mean, geom = "path") +
stat_summary(aes(color = Region), fun.data = mean_cl_boot,
geom = "errorbar", width = 0.1) +
stat_summary(aes(color = Region), fun = mean, geom = "point", size = 4) +
# geom_point(aes(color = Region), size = 2) +
theme_minimal() + ylab("cCDD") +
theme(text = element_text(size=12),
axis.text.x = element_text(angle=0))
ggplot(budbreak_cCDD3,
aes(x=Garden, y=cCDD, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path") +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 2) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Bud break in 2020 (cCDD)") + theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Edge"="steelblue","Margin"="green2"))
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
# diagnostic tests for the model
sjPlot::plot_model(budbreak_cCDD_mod1, type = "diag")
## [[1]]
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]
## [[2]]$Family
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]$Population
## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]$mBed
## `geom_smooth()` using formula 'y ~ x'
##
##
## [[3]]
##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'
tab_model(budbreak_cCDD_mod1)
| cCDD | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 80.26 | 80.17 – 80.35 | <0.001 |
| Garden [Maryland] | 13.54 | 13.37 – 13.71 | <0.001 |
| Garden [North_Carolina] | -14.78 | -14.91 – -14.64 | <0.001 |
| Region [Margin] | -0.01 | -0.13 – 0.11 | 0.875 |
| Region [Edge] | -0.00 | -0.10 – 0.09 | 0.951 |
|
Garden [Maryland] * Region [Margin] |
-0.06 | -0.35 – 0.23 | 0.698 |
|
Garden [North_Carolina] * Region [Margin] |
-0.01 | -0.20 – 0.17 | 0.888 |
|
Garden [Maryland] * Region [Edge] |
0.49 | 0.26 – 0.71 | <0.001 |
|
Garden [North_Carolina] * Region [Edge] |
0.42 | 0.27 – 0.57 | <0.001 |
| Random Effects | |||
| σ2 | 0.75 | ||
| τ00 Family | 0.00 | ||
| τ00 Population | 0.00 | ||
| τ00 mBed | 0.01 | ||
| τ11 Family.GardenMaryland | 0.57 | ||
| τ11 Family.GardenNorth_Carolina | 0.06 | ||
| ρ01 Family.GardenMaryland | |||
| ρ01 Family.GardenNorth_Carolina | |||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4757 | ||
| Marginal R2 / Conditional R2 | 0.994 / NA | ||